Coulomb Disorder Effects on ARPES and NQR Spectra in Cuprates 
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The role of Coulomb disorder, either of extrinsic origin or introduced by dopant ions in undoped 
and lightly-doped cuprates, is studied. We demonstrate that charged surface defects in an insulator 
lead to a Gaussian broadening of the Angle-Resolved Photoemisson Spectroscopy (ARPES) lines. 
The effect is due to the long-range nature of the Coulomb interaction. A tiny surface concentra- 
tion of defects about a fraction of one per cent is sufficient to explain the line broadening observed 
in Sr2Cu02Cl2, La2Cu04, and Ca2Cu02Ci2. Due to the Coulomb screening, the ARPES spec- 
tra evolve dramatically with doping, changing their shape from a broad Gaussian form to narrow 
Lorentzian ones. To understand the screening mechanism and the lineshape evolution in detail, 
we perform Hartree-Fock simulations with random positions of surface defects and dopant ions. 
To check validity of the model we calculate the Nuclear Quadrupole Resonance (NQR) lineshapes 
as a function of doping and reproduce the experimentally observed NQR spectra. Our study also 
indicates opening of a substantial Coulomb gap at the chemical potential. For a surface Cu02 layer 
the value of the gap is of the order of 10 meV while in the bulk it is reduced to the value about a 
few meV. 
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PACS numbers: 74.72.Dn, 79.60.Ht, 76.60.Gv, 73.20.At 



I. INTRODUCTION 

During last two decades, the Angle-Resolved Pho- 
toemission Spectroscopy (ARPES) has developed to be 
one of the most important methods in the physics of 
strongly correlated systems* 1 - Although the mechanism 
and physics behind the method is well understood, there 
are still issues remaining open to date. The large 
ARPES lincwidth observed in insulating parent com- 
pounds Si^CuC^C^, La 2 Cu04, and Ca2Cu02Cl2 is one 
of such problems. Already in the first experiments with 
Sr 2 Cu0 2 Cl 2 ^ it has been demonstrated that while the 
quasiparticlc dispersion is well described by the extended 
t-J model^£ the quasiparticlc width about 0.4 cV is 
difficult to reconcile with the predictions of this model 
alone. Similar broad spectra were later observed in 
La 2 CuOA!£ and in Ca 2 Cu0 2 Cl2^^. Interestingly, this 
broadening is not quite universal: while the linewidth 
in Sr 2 Cu0 2 Cl 2 and La 2 Cu0 4 is about 0.40-0.45 cV, it 
is only about 0.25 eV in Ca2Cu0 2 Cl 2 , possibly indicat- 
ing an extrinsic origin of the effect. Another important 
observation was made in Ref. 0: the quasiparticle lines 
display a Gaussian shape which is difficult to understand 
in terms of quasiparticle damping resulting typically in 
lineshapes of the Lorentzian form. 

Evolution of the ARPES lineshapes with doping has 
also been studied intensively i 6 ' 7 ' 8 ' 9 ' 10 At doping as small 
as 3%, the lineshape has already changed dramatically: 
a narrow peak of a Lorentzian shape is found to emerge 
from a shoulder-like broad background, with the intensity 
roughly proportional to doping. 

The doping dependence of the 63 Cu Nu- 
clear Quadrupole Resonance (NQR) spectrum in 
La2- a: Sr a ;Cu0 4 11 i 12 displays a totally different behavior. 
The NQR line is very narrow in the parent compound; 



the effect of doping is to broaden and shift the spectra 
to higher frequency. Since NQR is a local probe of hole 
density, the broad spectrum indicates a very inhomo- 
geneous hole density profile in the bulk of the sample, 
which is a consequence of the intrinsic disorder due to 
random La — * Sr substitutions 1 ^. 

An explanation of the broad ARPES lines in undoped 
cuprates was suggested in Ref. [t^. According to this sce- 
nario, the strong interactions between holes and optical 
phonons lead to the Franck-Condon broadening of the 
spectral functions. A detailed numerical study of a sin- 
gle hole in the t-J model coupled to optical phononsi 3 - has 
confirmed that by appropriate tuning of the holc-phonon 
coupling one can properly reproduce the ARPES spectra 
(see also Ref. [13 for a review). In the Franck-Condon 
broadening picture, multiple phonon subbands are gen- 
erated by a photoexcited hole and the observed broad 
ARPES line corresponds to the hole-phonon incoherent 
background. The narrow quasiparticle line still exists, 
but it is practically invisible due to strongly suppressed 
quasiparticle residue. 

While the electron-phonon mechanism enhanced by 
correlation effects in Mott system s 13 : 14 is able to explain 
broad ARPES lines in insulating cuprates, some ques- 
tions remain to be clarified. As noticed in Ref. ITU strong 
suppression of quasiparticle peak due to Franck-Condon 
mechanism implies also a drastic enhancement of a hole 
effective mass from its "bare value" m* w 2m e calculated 
within the t-J model. The resulting large mass polarons 
are then readily trapped by defects, e.g., by a negatively 
charged Sr dopant ion due to Coulomb attraction. For 
such a strong localization on atomic scales, the antifer- 
romagnetic order would survive up to a very high dop- 
ing level (similar to the case of Zn substitution), which 
contradicts the experimental data. In fact, the hole lo- 
calization length in a lightly doped La2Cu04 (x < 0.01) 
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is known to be about 10A corresponding to a moderate 
mass m* m 2m e (see Refs.[Hll3). Recent calculations^ 
suggest that nonlocal nature of electron-phonon interac- 
tion and longer-range hoppings may help to resolve the 
above difficulties. To reach a conclusive picture, however, 
further theoretical studies of the electron-phonon mech- 
anism in cuprates at finite density of holes are required. 

In the present paper we consider a different mechanism 
for broadening of the ARPES spectra. The mechanism 
is based on Coulomb disorder effects. There are two dis- 
tinct kinds of Coulomb defects under consideration. The 
first kind is related to the doping mechanism, where ran- 
dom La — > Sr substitutions create (negatively charged) 
Coulomb defects in the bulk. Bulk density of these de- 
fects is equal to the Sr concentration and hence equal to 
the doping level x. 

The second kind of defects are surface defects. We 
assume that cleaving the crystal creates some surface 
Coulomb defects that are unrelated to doping, and each 
surface defect has either positive or negative elemen- 
tary charge. The presence of surface defects (e.g., miss- 
ing surface ions) is physically plausible, and, in fact, 
they are observable with a scanning tunneling micro- 
scope (STM) 19 i 20 . Denoting density of positive defects 
by C, we assume that negative defects have the same 
density to fulfil the charge neutrality condition. While 
the value of C is material sensitive and not known a pri- 
ori, we will demonstrate that a concentration about a 
fraction of one per cent is already sufficient to explain ob- 
served ARPES broadening in undoped compounds. One 
may argue that the emergence of a narrow quasiparticle 
peak upon doping disfavors the disorder picture, since 
disorder is also enhanced upon doping.— However, in- 
teractions between holes induce nontrivial screening of 
impurities and dramatically reduce the effect of disorder 
on the ARPES spectra at finite doping. We find that 
the screening effects lead to the onset of narrow quasi- 
particle peaks already at doping as small as 1%. More- 
over, we will examine the effect of hole- hole interactions 
on the density of states (DOS), where we recover the 
well-known results of localization theor y 21 ' 22 and discuss 
their implication on transport properties of lightly doped 
La2- a: Sr a ;Cu04i 23 ' 24 i 25 Within the same model and ap- 
proximations, wc also address the disorder effects on 
NQR spectra and find a good description of the experi- 
mental data. 

Structure of the paper is as follows. In section II, we 
consider insulating undoped compounds and calculate ef- 
fect of surface Coulomb defects on ARPES spectra. In 
section III, wc consider doped compounds and calculate 
effect of bulk Coulomb defects on NQR spectra. Section 
IV highlights the effect of interactions on the bulk DOS. 
Evolution of ARPES spectra with doping is calculated 
in section V, where both surface and bulk Coulomb de- 
fects are taken into account simultaneously. In the case 
of doped compounds (Sections III, IV and V) one must 
consider screening of the long range Coulomb interaction 
by mobile holes, which is done numerically by perform- 



ing many-body Hartrce-Fock simulations. It is noticed 
that we do not consider a superconducting pairing in the 
present paper. Our primary goal here is to analyse role 
of the Coulomb disorder. Therefore, we concentrate on 
single particle properties and on the long range Coulomb 
interactions only. 

II. ARPES LINESHAPE IN AN INSULATOR: 
BROADENING BY SURFACE COULOMB 
DEFECTS 

We first consider the ARPES spectrum in the undoped 
insulating case, where a single hole is injected into the 
cleaved surface of, e.g., La2Cu04. Coulomb potential 
energy of the hole at position r due to interaction with 
surface defects is 

U(r)=J2 6g ' 2 ■ (1) 

, Wk - r i\ + a d 

Here e is charge of the hole (elementary charge) and qi = 
±e x q is charge of the defect. Eventually we will take 
q = 1, but now we keep q as a parameter for general 
Coulomb disorder. We assume that the defect is located 
at distance above the Cu02 plane, and r; is the 2D 
position of the defect. 

Note that we use electromagnetic units, I/47reo = f, 
and the effective surface dielectric constant is2£ 

e s = ^(e + l), (2) 

where e is effective bulk dielectric constant. According to 
Ref.0 the bulk dielectric constant is slightly anisotropic 
e c ~ 30, e a b ^40. In this case the effective bulk constant 
in Eq. © reads as^ e = ^/e a t,e c . As a representative 
value for cuprates, we will use e = 40 for the bulk dielec- 
tric constant throughout the paper. 

The lattice spacing of planar Cu's is set to be unity, 
ao ~ 3. 8 A — > 1, and hence the concentration of surface 
Coulomb defects C_ = C + — C is measured in units 
of the number of defects per Cu site. In the limit of low 
defect concentration, the potential energy constructed by 
Eq. {!]) varies slowly as a function of position r, and one 
can define a distribution function P(U) as the probability 
to find a given value of U. As C_ = C+, the average value 
of U is zero, U = J UP(U)dU = 0. It is instructive to 
calculate the root mean square deviation from zero, us 2 — 
U 2 — J U 2 P(U)dU . Squaring Eq. {T]) and averaging over 
random positions of defects we find 

u = ) q 2C / 2 , 3 

\e s a J J r 2 +a 2 d 

where L is the long distance (infrared) cutoff. Thus, 

u Q = V^4nC\n(j£j , (4) 

2 

V = — w q x IQOmeV , (5) 
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where V fixes the energy scale. For purely random distri- 
bution of Coulomb defects the infrared cutoff L is equal 
to the radius of the incident photon beam (assuming that 
the radius is smaller than the sample size). In principle 
one can also imagine that due to some reason related 
to the cleaving process there is a long distance correla- 
tion between Coulomb defects. In this case the infrared 
cutoff is equal to the corresponding correlation length. 
Fortunately the dependence of ujq on the infrared cutoff 
is very weak. For instance, loq for L = 1mm differs from 
that for L = \\im only by 30%. It is worth mentioning 
that the infrared logarithmic divergence of the integral 
in Eq. ([3]) is a consequence of the long range nature of 
the unscreened Coulomb interaction. 

It is also instructive to study the spatial correlator of 
the potential, (U(r)U(Q)). Clearly, there is a structure in 
the correlator at a distance about average separation be- 
tween defects, r ~ 1/y/C. However, the most interesting 
behavior is at distances 1/VC <C r <C L. A straightfor- 
ward calculation similar to ([3]) yields that in this regime 



(U(r)U(O)) 
(U*(0)) 



In I 



In 



(6) 



thus the potential varies slowly at the scale comparable 
with the infrared cutoff. 

Now we calculate the entire distribution function of the 
potential, 



P(u) = (6{w - U(r))) 



(7) 



where (...) denotes averaging over the observation point 
or, alternatively, averaging over distribution of defects. 
We choose to put the observation point at the origin 
and perform averaging over distribution of defects, hence 
Eq. (P) becomes U(0) = J2i W ~ J2j JT> wnere * enu ~ 
merates positive defects and j enumerates negative ones, 
and R = \/ r 2 + a 2 , is the distance measured in units of 
lattice spacing. From Eq. (J7J) one obtains 



PW) = s 



+ 



dte 1 



-E-t TT +ijht\ 



(8) 



All defects are distributed independently, therefore 

(Qe-^n^) - (e-^We***)?- , (9) 



where N + and N_ are total numbers of positive and neg- 
ative defects and (...) r denotes averaging over the defect 
position. Let us denote by N the total number of sites 
in the square lattice, therefore 



where 



y 



I±= j dfr ( 1 - cxp { Ti—t 



(10) 



(11) 



When deriving (fTO]) we keep in mind that J cPr = 1. 
Hence 



I± 



N± 



-CI± 



(12) 



Here we have taken into account that concentration of 
defects is C = N+/N = N-/N. Hence, Eq. © is trans- 
formed to 



e + ^*) = e -c(i + +i-) 



cxp <^ ~2C d 2 r\l- cos 



Vt 
~R 



(13) 



To evaluate the integral with logarithmic accuracy, we 
expand cos at <C 1 , obtaining 

r-L 



2C I d 2 r | 1 - cos 
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(14) 



Substituting Eqs. (jUI) and (JTHJ) into Eq. (Hj) and perform- 
ing integration over t, we find the Gaussian distribution 
for the potential: 



P(w) = 



2ttujo 



(15) 



where luq is given by Eq. ^j. Hence the half width of the 
potential distribution is 



T = v/8 In {2)uj = 4j27rln(2)ln 



a,! 



xVCxV . (16) 



A numerical simulation, that statistically includes 100 
defect configurations for system of size — > 10, shows a 
remarkable consistency between the potential generated 
by Eq. (Q} and its analytical distribution, Eq. (fT5|) . with 
the width (fTB")) . Note that one should be cautious about 
the numerical value of ad- Since it is the effective short 
range cutoff of the Coulomb potential, it must also in- 
clude the size of Zhang-Rice singlet which is about one 
lattice constant. The precise value of ad will be discussed 
in the Sec. Ill, but now we take = 3.8A — ► 1. The 
width of potential distribution in Sr2Cu02Ci2, La2Cu04, 
and Ca2CuC>2Cl2 can then be estimated by choosing 
g=l and L = Imm, which results in 



r 
r 



250 meV 
450 meV 



at 
at 



C 

c 



0.002 
0.006 



(17) 



We suggest that Eq. ((TSJ) is exactly the ARPES 
line broadening function with the half-width given by 
Eqs. (fT6|) and (jTTJ) . This can be understood as follows: 
In the photoemission process a single hole is injected in 
the top layer of the insulator, and there are two mecha- 
nisms for the line broadening in the presence of disorder. 
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The first one is the direct scattering of the hole from in- 
dividual defect, which is the short range mechanism and 
therefore its contribution to the broadening is propor- 
tional to the first power of concentration of defects. We 
expect that this mechanism is negligible because of the 
low concentration of defects. The second mechanism is 
due to the fact that different holes arc injected in differ- 
ent parts of the sample which have different potentials. 
Spectral functions are then broadened due to the poten- 
tial distribution ([15)) . Contribution of this mechanism 
to broadening is proportional to the square root of con- 
centration of defects, and, moreover, it is logarithmically 
enhanced, see Eq. (fl"6|) . 

The above picture is supported by numerical calcu- 
lations that we going to discuss now. To construct a 
model that can properly describe a hole motion in un- 
doped cuprates we first notice that there are different 
length scales in the problem: (i) The scale of the order 
of 1-2 lattice spacing. Strong correlations, such as exci- 
tations of multiple virtual magnons, occur at this scale. 
(m) The scale about average separation between Coulomb 
defects ~ 1/y/C. Scattering from defects takes place at 
this scale. (Hi) The scale of 1/yC « r < L, where loga- 
rithmically enhanced variations of the potential develop. 
Regarding the point (i), we do not treat here the strong 
correlations explicitly, but adopt the effective dressed 
hole dispersion after quantum fluctuations at short dis- 
tances (i) are included. It is known that dispersion of the 
dressed hole has minima at points (±7r/2, ±7r/2), and is 
approximately isotropic around these points^. The band 
width of the dressed hole is about 2J, where J m 130 
meV is the superexchange in the t-J model, although we 
do not directly employ the t-J formalism. Hereafter we 
set energy units 



J = 130 meV -> 1 . 



(18) 



To imitate dispersion of the dressed hole we consider spin- 
less fermions on a 2D square lattice. The Hamiltonian 
reads as 



H, 



(19) 



where c| is the hole operator at site i, and t" denotes 
the next-next-nearest-neighbor hopping on the square 
lattice. The Hamiltonian (fTS)) yields the following dis- 
persion 



£k = 2t" (cos2k x + cos2k y 



(20) 



The dispersion is isotropic around minima at points 
(±7r/2, ±7r/2) as shown in Fig.Q] We choose t" = 0.25 to 
reproduce the realistic hole bandwidth as obtained from 
the t-J model. Notice that in the original t-J model for- 
malism, there are four half-pockets inside magnetic Bril- 
louin zone, and each pocket has two pseudospins^ in the 
present model we consider four full pockets inside the full 
Brillouin zone with spinlcss fermions, hence the number 
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FIG. 1: Dispersion minima of the spinless fermion generated 
by Hamiltonian (11911 , 



of charge degrees of freedom is exactly the same. Clearly 
our model does not have a momentum dependence of the 
quasiparticle residue, especially its suppression outside 
of magnetic Brillouin zone due to strong correlations £ 
However, the important point is that the Coulomb in- 
teraction remains unchanged whatever the value of the 
residue is. This is because the charge is conserved even 
though holes are heavily dressed. 

The hole-defect interaction due to Coulomb potential 
in Eq. fTJ) reads, after we set J — * 1 and ao — ► 1, as 
follows: 



Hfi-d — 



E 



Qi 



with a dimensionless "charge" value 
Qi = ±j « ±1.5 



c\ci , 



(21) 



(22) 



The superscript "s" stands for "surface" . This yields the 
full Hamiltonian 



H — H t + Hh-c, 



(23) 



which can be easily diagonalized on a finite size cluster 
where positive and negative defects with concentration 
C each are randomly distributed. The ARPES experi- 
ments measure the electron spectral function which can 
be calculated exactly using the cluster eigenstates and 
eigenenergies. Denoting the hole energy as e and the 
electron energy as u>, we have w = — e, and hence 



A(k,e) = £|<k|n)| a *(e-.E n ) 



TV 



(e-E n y+r] 2 



A(k,w) = A(k,-e) , 



(24) 



where a n (ri) = (i\n) is the coordinate representation of 
the n— th cigenstate, H\n) = E n \n), and rj = 0.01 is the 
artificial broadening of discrete energy spectrum. Note 
that the artificial broadening rj = 0.01 = 1.3 meV is much 
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smaller than any physical contribution to the width con- 
sidered in the present work. Therefore r\ is just a techni- 
cal mean to stabilize the numerical procedure and does 
not contribute to physical broadening. We perform diag- 
onalization in a 36 x 36 cluster with periodic boundary 
conditions, where the distance |R; — r^| is chosen to be 
the shortest distance on the torus. A statistical averaging 
over 100 disorder configurations is performed. 



the radius of the incident photon beam. In this case fo- 
cusing/defocusing of the beam changes the width accord- 
ingly. For instance, the radius variation from L = l//m 
to L = 1mm gives rise to a sizable variation of the width 
by 30%. It can be noticed, however, that the width is 
most sensitive to the amount and the nature of defects, 
so systematic studies of the width variation as a function 
of the surface quality might be useful. 
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FIG. 2: (Color online) ARPES spectral function (J24J) for 36 x 
36 cluster averaged over 100 random realizations of Coulomb 
disorder. Concentration of defects is C = 0.6%. Values of 
momenta are k = (n /2,7r/2), (0,7r/2), and (0,0). Smooth 
lines correspond to energy dispersion u — —e^ [with tk given 
by Eq. ([20])], broadened according to Eqs. 1(15). (Tl6)) , and (0). 
In these equations, we set L = and a<j = 1. 

Fig. [2] shows the resulting ARPES spectral function 
at three different momenta k = (tt/2, 7r/2), (0,7r/2), 
(0, 0) together with the analytical distribution (fT5|) of 
the Coulomb potential P(ui). In this pedagogical exam- 
ple where the cluster size is certainly smaller than size 
of experimental samples, we choose C = 0.6%. There 
is an excellent agreement of the line shape and the line 
width between P(uj) and A(h, uj) in all three momenta. 
This confirms our statement that Eq. (fT"5|) describes a 
Gaussian broadening of ARPES spectra in an insulator 
due to surface Coulomb defects, and implies that accord- 
ing to Eq. (|17|). the concentration of 0.6% surface defects 
(for one charge species) is sufficient to explain the ob- 
served broadening in La2CuC>4, and 0.2% is sufficient for 
Ca2Cu02Cl2. It is worth mentioning that this broaden- 
ing is a fairly general mechanism that works not only for 
Mott insulators, it is also valid for usual band insulators. 

How one can check experimentally validity of the sug- 
gested broadening mechanism? There are a few "han- 
dles" in Eq. (fT6|) : the infrared cutoff L, the ultraviolet 
cutoff ad, the density C and the potential V of the de- 
fects. One can vary a<j by probing deeper layers (e.g., 
using different photon energies): for the top layer Od = 1, 
for the second one ad = 4.4, for the third layer ad = 7.8, 
etc. Therefore the width for the second layer must be by 
5-6% smaller than the width for the top one, the width 
for the third layer is further reduced by 2-3%. These are 
rather small effects. Another possible "handle" is the in- 
frared cutoff L. It is natural to assume that L is equal to 



III. SPECTRUM OF COPPER NUCLEAR 
QUADRUPOLE RESONANCE 

In this section we examine the effect of Coulomb dis- 
order on the NQR spectrum. Doping dependence of 
63 Cu NQR in La2-- L Sr, E Cu04 has been studied in de- 
tail in Refs. [HI and [H, which show that in the undoped 
parent compound the NQR spectrum is comprised of a 
very narrow line (centered at frequency 33.05 MHz at 
T=600 K). The spectrum is shifted to higher frequency 
upon hole doping, with a line width roughly proportional 
to doping. Since NQR is a local probe of real space hole 
distribution,— the broad spectrum indicates a very inho- 
mogencous profile of hole density^ Another interesting 
feature is that NQR spectra in doped samples show a 
double structure: a secondary hump (the "B-line") ap- 
pears at a frequency higher than the broad main line. 
The origin of the B-line is attributed to the Cu sites 
that are directly underneath the Sr substitutions. We 
will show the experimental data and compare it with our 
results later in this section. 

The NQR spectrum is obtained by calculating the hole 
density distribution which is spatially nonuniform due 
to disorder. The model used in the previous section is 
modified here as follows. Since NQR is a bulk sensitive 
measurement the surface defects are not relevant, and 
disorder effect is solely due to the Coulomb potential of 
randomly distributed Sr-dopants. The total number of 
out-of-planc Sr ions is equal to that of the holes, as sug- 
gested by the doping mechanism of La2-^Sr x Cu04, and 
each Sr-defect brings about a negative charge. The con- 
centration of negative defects is therefore equal to doping, 
and no positive defects are present. Further, the strength 
of Coulomb interaction is reduced comparing to the sur- 
face case, because the bulk dielectric constant is larger: 
e ~ 2e s , see Eq. <j2j> . We denote the effective dimension- 
less charge in bulk as 

Q « ~ 0.75 , (25) 

which is half of the surface dimcnsionlcss charge Q s , 
Eq. |22|) . The interaction of a hole with Sr ions is then 

l,i 

U u = Q - . (26) 
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Similarly, the Coulomb interaction between holes is de- 
scribed by 



y/\n-rj\ 2 + ai ' 



(27) 



where we use the same cutoff ad = 1 to represent the size 
of the Zhang- Rice singlet. 

We further consider the effect of multilayer screening 
in bulk of I^-xSr^CuO^ that contains a periodic struc- 
ture of CuC>2 layers along c-axes. The electric field of 
a charge in a particular layer is substantially screened 
and deformed by the other layers, as shown schematically 
in Fig. [3l due to their large polarizability. In principle, 
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-Q 



E 



Ui 
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\ V|R-*-ri| 2 +"(2W, 
1 



Ti -r 3 \ 2 + a' d 

2(-l) n 
^ ^\v t -vtf + {2ndY 



(28) 



where d = 13. 2A — > 3.5 is the separation between lay- 
ers. The crossover between these two limiting cases corre- 
sponds to the situation when the in-plane dielectric con- 
stant due to holes is equal to the ionic dielectric constant, 
which takes place at x ?« Since we are interested in 

the range of x > 1%, relevant to the NQR data discussed 
below, Eq. (f2"B"]) is applied to all finite doping cases in our 
simulations. Accuracy of this approximation is some- 
what questionable at x = 1%, but it is fairly reasonable 
at higher dopings. 



0.2r 



0.15- 



0.1 



a image charge 

FIG. 3: Screening of in-plane Coulomb interaction by other 
layers. Dashed lines shows electric field of an in-plane charge 
bent due to a large polarizability of "other layers" . 



one needs to perform a self-consistent calculation that 
includes multiple layers to account for this effect. Un- 
fortunately, such a calculation is too expensive compu- 
tationally. However, one can consider the following two 
limiting cases where analytical descriptions are available. 
The first case is that at extremely small doping, the po- 
larizability of other layers is negligible, hence we recover 
the single layer formulae in Eqs. (|2"rH) and (f2"T|) . The sec- 
ond case is that at sufficiently large doping, the high po- 
larizability implies that the electric field generated from 
one layer is practically perpendicular to the surface of 
its nearest layers as it is shown in Fig. [3] In this case, 
one can simply apply the method of image to account 
for the screening due to other layers, by assuming that 
each layer is placed between two highly polarizable me- 
dia. The interactions described in Eqs. (|2"rH) and (|2~T|) are 



05 1 1 1 1 1 1 1 1 

u 0.6 0.8 1 1.2 1.4 

a d 

FIG. 4: The hole-Sr binding energy as a function of the short- 
range cutoff ad- We recall that we set J = 1. The interlayer 
screening is taken into account according to Eq. (I28|l . 

Here we give more details about the choice of a<j, under 
the condition that the multilayer screening has been ac- 
counted for via Eq. (|28p . The main reason for this short 
range cutoff is the size of the Zhang- Rice singlet, which 
a priori yields ad ~ 1. On the other hand, the value of 
ad should properly restore the binding energy of a sin- 
gle hole trapped around a Sr ion, which is known to be 
£b ~ 10 mcV <~ O.U." 9 To estimate e&, we assume that 
there exists a large enough doping range where multi- 
layer screening takes place via Eq. (|28|) , while the doping 
is still small enough that the in-plane hole-hole inter- 
action can be ignored, and diagonalize the Hamiltonian 
H = Hi + Hh-Sr with only one Sr present. The result- 
ing binding energy versus ad is shown in Fig. 01 where 
we found that ad = 1 indeed gives the correct binding 
energy. For extremely low doping x <C 1%, we adopt the 
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unscreened potential Eq. (|26|) without considering other 
layers, and found that = 1 gives e b ss 0.23J « 30 
mcV. This demonstrates the importance of the multi- 
layer screening at doping x > 1%. The value = 1 is 
adopted throughout this work. 
The full Hamiltonian 



H = H t + H, 



h-Sr 



Hi, 



(29) 



can be diagonalizcd by the following analysis. The di- 
mensionless parameter that characterizes the strength of 
interaction in a 2D Coulomb gas is^ 



0.36 



eh' 



(30) 



We see that even at x = 2% the value of r s is still small 
r s f» 1.4. Moreover, the multilayer screening introduced 
in Eq. (|2"8j) further reduces this value to r s — > 1. There- 
fore, we are safely in the weak coupling regime where the 
Hartree-Fock treatment is adequate. Notice that the hole 
dynamics are certainly strongly correlated at the length 
scale about a few lattice spacing. These are Hubbard or 
t-J model correlations which result in the dispersion of 
dressed holes, Eq. ([20]) . Here the term "weak coupling" 
refers to the long range Coulomb interaction between 
holes at the length scale > where the effect of 

the short-range strong correlations is already taken care 
of by adopting the dispersion (j2"0)) . The Hartree-Fock de- 
composition is then applied to the hole-hole interaction 



H,, 



u, 



(31) 



The diagonalization is again done in a 36 x 36 cluster 
with periodic boundary conditions (torus). Expectation 
values (c\ci) and (c\cj) can be calculated by 

(4a) = £K(i)| 2 /(£n) , 

n 

<ck> = E Q « W* a ™ U)f(E n ) , (32) 



where 



f(E n ) 



(33) 



is the Fermi-Dirac distribution. 

To determine the macroscopic chemical potential in 
our simulation, we adopt the following procedure. The 
chemical potential at either zero or finite temperature 
in each defect configuration is determined via the charge 
neutrality condition, i.e. the number of negative defects 
is equal to the number of holes. The average value of the 
chemical potential, denoted by [i, is then calculated out of 
100 defect configurations taken. We then shift the entire 
energy spectrum of each particular configuration in such 
a way that the chemical potential of the configuration 



x=2% 



T=0K 




FIG. 5: (Color online) Plots of the hole density in a CuCh 
layer deep in the bulk, for a particular realization of random 
Sr positions, at doping x = 0.02 and two different tempera- 
tures T = and 600 K. The length scale 5 nm (~ 13ao) is 
shown. 



is equal to this mean value ^, which is the macroscopic 
chemical potential. 

Hole density plots m = (eje,-) for particular realiza- 
tions at x = 0.02 and x — 0.07 are shown in Figs. [5] 
and [51 respectively, for two different temperatures T = 
and T — 600 K. One sees very inhomogeneous density 
profiles, with a characteristic length scale of the order 
of a few nanometers. Similar nanoscale charge inhomo- 
geneities, that closely resemble the STM images of un- 
derdoped cuprates (see, e.g., Refs.[l9l.[3ll). have been also 
reported in previous studie a 32 i 33 ' 34 . Interestingly, we find 
that increasing temperature substantially reduces the in- 
homogeneity (compare upper and lower panels in Figs. [S] 
and [5]) ■ This has consequences for NQR spectra which 
we address in the following. 

The NQR frequency at a particular site i is related to 
the hole density m by2£ 



(33.05 + 19 • m)MHz 



(34) 



Thus, the entire NQR spectrum, which effectively sums 
over all sites in the sample, is proportional to the prob- 
ability distribution of the hole density V(n), up to a 
constant shift 33.05 MHz corresponding to the parent 
compound^. Calculated NQR spectra for T = and 
T = 600 K are presented in the upper and in the middle 
panels of Fig. [71 respectively. The experimental plots of 
Refs. UH and \lU for I^-xSr^CuC^ at different dopings 
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x=7% 



T=OK 




32 



34 



36 



38 



40 



FIG. 6: (Color online) Plots of the hole density in a C11O2 
layer deep in the bulk, for a particular realization of random 
Sr positions, at doping x = 0.07 and two different tempera- 
tures T — and 600 K. 



and T=600 K arc shown in the lower panel for compar- 
ison. One sees that the shift and the broadening of the 
spectrum upon doping are well reproduced by the theory, 
with the lincwidth consistent with experimental data. 
Surprisingly, the agreement is reasonable even at 16% 
doping which is certainly too large for our low doping 
theory. The mechanism that narrows the linewidth with 
temperature^ is also understood: increasing tempera- 
ture reduces the spatial inhomogencity of density profile, 
as shown in Figs. O and [51 which results in a narrower 
probability distribution V{n), and hence the narrower 
NQR spectrum (compare the upper and the middle pan- 
els in Fig. [7]) . Although our simple model does not take 
into account the direct action of the Sr ion Coulomb field 
on the nearest Cu nuclei, which is believed to be the main 
origin of the high frequency " B-line" ^ in our numerics 
we do see a shoulder-like structure emerging at high fre- 
quency. A detailed investigation shows that this struc- 
ture is associated with holes that are trapped around 
local potential minimum due to occasionally clustered 2 
to 3 Sr defects. We suspect that these trapped holes, al- 
beit not the main reason for the " B-structure" , can have 
a certain contribution to it. The most important point 
here is that the main line is well understood in our model, 
which properly captures the hole density distribution and 
the screening effects in presence of Sr defects. 



CZ 
Z3 



CO 



f0.5 

£= 



2% 4% 7% 11.5% 16% 



T=0K 




Z3 



CO 



f 0-5 
c 

CD 



2%4%7% 1 .5% 16% 



T=600K 



0% 




2%4%7%11.5% 16% 



T=600K 




36 
f(MHz) 



FIG. 7: (Color online) The calculated (upper and middle pan- 
els) and experimental (lower panel) 63 Cu NQR spectra in 
La2^ a; Sr a ;Cu04 at different doping levels x= 0.02, 0.04, 0.07, 
0.115, and 0.16. The narrow line at / = 33.05MH^ corre- 
sponds to the parent compound, x — 0. Apart from broad- 
ening of the spectra, doping results also in a high-frequency 
structure similar to the "B-line" observed in the experiment. 
Note that linewidth at T=600 K is smaller than that at zero 
temperature. 



IV. DENSITY OF STATES AND ANDERSON 
LOCALIZATION IN Cu0 2 LAYER IN THE BULK 

In this section we address the issue of the bulk DOS 
in presence of disorder, and relation of the DOS to the 
normal state dc electrical conductivity. It is known that 
La2_ 2; Sr 2 ;Cu04 exhibits the variable range hopping con- 
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ductance at small doping x < 0.055j 17 ' 23 . This indicates 
a strong localization of holes by Sr Coulomb potential. It 
has been suggested that the onset of superconductivity 
at x > 0.055 is due to percolation of the bound states^ 
The hole density plots shown in Figs. [5] and [6] support this 
suggestion: the hole density at x = 0.02 (T = 0) vanishes 
in large areas of the system, signaturing a highly local- 
ized density profile, while the density at x = 0.07 (T = 0) 
is nonzero practically everywhere in the system, so wave 
functions are highly overlapped. 

To study the problem in more detail, we have calcu- 
lated the 2D density of states, p(e), via the standard 
definition: 

p(e) = ^X>(e-£„). ( 35 ) 

n 

Fig. [8] shows the DOS calculated using the eigenencrgics 
of Eq. (|29|) and fixing the chemical potential as described 
in the previous section. One sees clearly a full reduc- 
tion of DOS at the chemical potential as expected by the 
Coulomb gap theory in 2D systems^. The size of the 
gap is about Ac ~ 2.0 — 2.5 meV at x = 0.02 and it 
decreases to about 1.5 — 2.0 meV at x = 0.07. We de- 
fine Ac such that the total width of the gap structure in 
DOS is 2Ac as indicated in the lower panel of Fig. [10] 
Importantly, the gap smoothly evolves through the per- 
colation point x = 0.055. This implies that for single 
particle dynamics the system remains an Anderson insu- 
lator even after percolation. Probably at even larger dop- 
ing, x > 0.1, the Coulomb gap smoothly evolves to the 
logarithmic reduction of DOS corresponding to the weak 
localization theory^ Unfortunately, we are not able to 
trace this crossover because the relatively small size of 
the cluster, 36x36, limits our accuracy of the gap cal- 
culation at the level ~ 1 meV. The insulating behavior 
across the percolation point x = 0.055 obtained in the 
present calculation agrees perfectly with the experimen- 
tal dat a 24 i 25 where the in-plane resistivity, measured in 
a very strong magnetic field that destroys superconduc- 
tivity, shows an insulating behaviour below ~ 50 K for a 
wide doping range up to x « 15%. 

The DOS displayed in Fig.[5]exhibits oscillations above 
chemical potential. These oscillations is a byproduct of 
the finite size of the cluster. Maxima of the DOS cor- 
respond to degenerate states with dispersion (|20|) on the 
36 x 36 torus. The oscillations must certainly disappear 
in the thermodynamics limit. However, oscillations of 
this kind also have an interesting physical meaning. In 
particular, the x = 0.02 plot in Fig. [8] indicates that 
while the quantum states near the chemical potential are 
strongly localized, the high-energy states well above the 
chemical potential are quite extended with a mean free 
path exceeding the size of the cluster used. Similarly, 
the smeared oscillations in Fig. [5] x — 0.07 indicate that 
the mean free path is comparable with or less than the 
cluster size. 
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FIG. 8: The hole DOS in Cu0 2 layer deep in the bulk at 
dopings x — 0.02 and x = 0.07. DOS vanishes at the chem- 
ical potential indicating the Coulomb gap which gradually 
decreases with doping. 



V. EVOLUTION OF ARPES SPECTRA WITH 
DOPING, AND DENSITY OF STATES IN THE 
SURFACE Cu0 2 LAYER 

To study the evolution of ARPES spectrum upon dop- 
ing, we apply the above Hartree-Fock treatment and in- 
tcrlayer screening picture to the top Cu02 plane near the 
surface. The ARPES lincwidth is determined now by a 
combined action of two types of disorder: the surface 
defects, as described by Hh~d in Sec. II, and a ran- 
domly distributed Sr-dopant ions described by Hh-Sr 
in Sec. III. To be specific, we fix here the concentra- 
tion of positively/negatively charged surface defects to 
be C+ = C_ = C = 0.6%, as required to fit the ARPES 
linewidth in La2Cu04, see Eq. (|17l) . We assume that 
concentration C is independent on Sr-doping since it is 
determined by the surface properties unrelated to dop- 
ing. Total concentration for negatively charged defects is 
then x + C = .x + 0.6%, counting both Sr-dopants and the 
negatively charged surface defects. Dimcnsionlcss charge 
Q s is again described by Eq. (f2"2"|) , which is twice of its 
bulk value Q due to the reduced dielectric constant ([2]) 
on the surface. Consequently, this enhances the disorder 
effects on the surface as we will sec below. The multilayer 
screening of the interactions on a cleaved surface is also 
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different from that in the bulk, Eq. ([28]) . because even 
though other planes are still considered as highly polar- 
izablc, and hence the method of image is still valid for 
x > 1%, the cleaved surface is now considered as located 
at a distance d above a polarizable media, see Fig. [9l in- 
stead of being sandwiched between two polarizable slabs. 
Collecting all these effects, we have 



H, 



h-d 



H, 



h-Sr 



£Z4 



c]ci 



(36) 



where a disorder potential, originating either from sur- 
face or Sr defects located at position R/, is given by 



±Q S 



1 



-1 

v/|R,-r. t | 2 + (W 



(37) 



The difference between Hh-d and Hh-Sr is only in the 
sign of Q s : Hh-Sr potential is always attractive and has 
charge —Q s , while its sign in Hh-d depends on charge of 
the surface defect which can be cither positive or nega- 
tive. Similarly, the hole-hole interaction reads as 



Hint ^ ^tlijCjCjCjCj ) 



(38) 



Ui 



cr 



i 



V\^~r^ + ai VN-r j | 2 + (2d)V 

We then diagonalizc the full Hamiltonian 

H = H t + H h - d + H h sr + H in t (39) 

following the Hartree-Fock treatment described in 
Eqs. ([31) to ([33]) . The procedure described in Sec. Ill 
is again applied to determine the macroscopic chemical 
potential. 

top "Hartree-Fock" layer 

» i other layer 



which is due to enhanced dimcnsionless charge Q s : a 2D 
Coulomb gap scales approximately as Q 2 , see Ref. l2ll . 

There are no sizable oscillations of the DOS in Fig.fTUl 
This means that the hole mean free path even well above 
the chemical potential is much smaller than the cluster 
size. 

Notice that values of the chemical potential on the sur- 
face, Fig. [T01 are somewhat different from those in the 
bulk, Fig. [5] This difference is a byproduct of our ap- 
proximations. A tiny surface charging energy and/or a 
tiny surface lattice deformation due to La — > Sr substi- 
tutions can tune up the surface chemical potential from 
its bulk value. Due to these effects, which arc not taken 
into account in the present model, we cannot compare the 
calculated chemical potential with experimental values. 
Although these effects can shift the chemical potential 
and overall energy scales, they do not influence the wave 
functions and the shape of DOS. 




x=7% 




m image charge 

FIG. 9: Screening of in-plane Coulomb interaction in the top 
layer by other layers underneath. Dashed lines shows electric 
field of an in-plane charge bent due to polarization of "other 
layers" . 

First, we discuss DOS obtained from the numerical 
diagonalization. Fig. [TU] shows the calculated DOS for 
the top CUO2 layer, where one clearly sees a Coulomb 
gap of the order of Ac ~ 0.1J ~ 10 meV opening at 
the chemical potential. Comparing the result with the 
bulk DOS shown in Fig. [51 we see a significant difference 



FIG. 10: The hole DOS in the surface CUO2 layer at dopings 
x — 0.02 and x = 0.07. DOS vanishes at the chemical po- 
tential indicating the Coulomb gap Ac- The size of the gap, 
Ac ~ 10 meV, is larger than that in the bulk (see Fig. [8}. 



Now, we turn to the ARPES spectra at finite dopings. 
The spectral functions are again calculated by using the 
eigenstates and eigenenergies given by exact diagonal- 
ization. Notice that at finite doping and zero tempera- 
ture only states above chemical potential are considered. 
This is because we are working in the hole representa- 
tion, while ARPES spectrum is associated with the elec- 
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tron spectral function, therefore only states occupied by 
electrons E n > fj, should be summed over in Eq. It 
is also convenient to shift origin of the momentum to the 
dispersion minimum (tt/2, 7r/2), 



(7r/2,7r/2) 



(40) 



and present the spectral function in terms of p. In our 
model A(jp,lo) is roughly symmetric around the disper- 
sion minimum, A(p,u>) ~ A(|p|,u;). In addition, because 
of the 36 x 36 finite cluster size, we can only calculate cer- 
tain discrete values of momentum, p = yg(m, n), where 
m,n are integers. On the other hand, the doping level x 
can be varied continuously. In Fig. [Til we present spec- 
tral functions calculated for the following momenta in the 
nodal direction 



Po 
Pi 

P2 
P.! 



= o , 



18' 18. 
2tt 2tt\ 

18' 18 ) 



37T 3tt 
18' 18 



(41) 



and for x = 0.01 — 0.11. Surprisingly, we see very nar- 
row lines with the width of the order V ~ 0.2 J ~ 30 
mcV, in spite of the very strong disorder. This is cer- 
tainly due to the Coulomb screening of both the surface 
defects and Sr dopant potentials. This residual "small" 
width r ~ 30 meV is due to the hole scattering from the 
residual (locally unscreened) part of the random poten- 
tial. Remarkably, the residual width is quite universal, 
it is practically independent of doping. On a qualitative 
level, this (somewhat unexpected) observation can be un- 
derstood as a result of two opposite trends: On the one 
hand, more Sr-doping increases amount of disorder but, 
at the same time, adding more holes reduces an effective 
scattering amplitude on each Sr-defect because of better 
screening. Within the present low energy effective theory, 
which is valid up to energies E ~ J ~ 100 — 200meV, the 
residual width is also energy/momentum independent. 

Finally, the evolution of the small Fermi surface upon 
doping can be considered. We first notice that in the 
homogeneous case the dispersion is given by Eq. ([20]) . 
Hence the Fermi momentum pp and the doping x are 
related as 



(42) 



if doping is small and the dispersion is roughly parabolic. 
Therefore, pp — pi at x w 0.02 and pp = Pi at x w 0.08. 
Interestingly, the relation (|42p remains qualitatively cor- 
rect even in the presence of strong disorder. This can 
be seen by the following analysis that extracts the Fermi 
momentum: We know that without disorder, the ARPES 
intensity at p below pp vanishes at any finite doping, 
since the momentum is inside hole Fermi surface and no 
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FIG. 11: (Color online) The electron spectral function A(\a, u) 
at doping levels x = 0.01, 0.03, 0.05, 0.07, 0.09, 0.11 calculated 
for momenta po, pi, P2, P3 specified in Eq. (|4ip . The 
Coulomb gap Ac is highlighted in the x = 0.11 plot. 



electron with this momentum can be excited. Therefore, 
the ARPES intensities at this momentum (red curves) 
in Fig. [TJJ are nonzero only because of disorder. The 
intensity decays very quickly when doping is increasing. 
We found that maxima of spectral functions are never 
exactly at the chemical potential, as one would expect 
in a system without disorder. Instead, we see that the 
maximum of each line gets closer to the chemical poten- 
tial as doping is increased, but stops at an energy scale 
Ac ~ 0.1J ~ 10 meV below the chemical potential. This 
is clearly due to the Coulomb gap opening in the DOS, as 
shown in Fig. 1 101 which suppresses the spectral function 
within /i — Ac < e < p and shifts the maximum. The 
maximum of the p = p\ line (green curves) approaches its 
rightmost position p — Ac at x w 0.03, so at this doping 
we say the " Fermi surface" crosses the momentum p = pi. 
Similarly, the maximum of the p = P2 line (blue curves) 
approaches its rightmost position at x ks 0.09, hence the 
Fermi surface crosses p = P2 at this doping. Following 
this procedure, we can identify the "Fermi momentum" 
at each doping, which can be compared with the experi- 
ments, although only discrete values at p = jg{m, n) can 
be identified. 

It should be stressed that our aim here is to study the 
Coulomb disorder and its influence on the quasiparticlc 
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peak. Strong (Hubbard or t-J model) correlations are 
taken into account only via effective dispersion (f2"0)) of 
the holon. In all other respects we disregard the spin 
degrees of freedom. Therefore, we do not reproduce the 
asymmetry of the ARPES spectral function A(p, ui) with 
respect to the boundary of the magnetic BriUouin zone, 
and also strongly underestimate the ARPES intensity be- 
low the quasiparticle peaks, seen in the experiment as a 
pronounced "hump" structure. In spite of these draw- 
backs, the theory allows us to address the issue of evolu- 
tion of the quasiparticle peak with disorder/doping. In 
particular, our theory explains why the ARPES lines in 
doped cuprates are relatively narrow in spite of the very 
strong Coulomb disorder.— iSuLitt 

Another interesting result obtained in this section is 
the predicted Coulomb gap ~ 10 mcV in density of states 
that could be observed by surface sensitive probes like the 
STM and ARPES. In fact, the gap features of this scale 
are present in the STM data for underdoped cuprates 
(see, e.g., Ref. Hl|), and are typically attributed to the 
(local) pairing gap; the Coulomb gap might be an addi- 
tional origin of these low-energy structures. 

VI. CONCLUSIONS 

In this paper, a comprehensive study of Coulomb 
disorder effects in undoped and lightly-doped cuprates 
is performed, and the main results can be summarized 
as follows. 

1. We have demonstrated that a very small amount of 
surface Coulomb defects leads to a dramatic broadening 
of ARPES spectrum in insulators. In particular, a 
concentration of defects about just a fraction of 1% is 
sufficient to explain observed ARPES line widths in 
La2Cu04 and Ca2Cu02Cl2- The broadened spectrum 
displays a Gaussian shape, consistent with experiments, 9 
In the end of section II we have discussed possible 
ways to check the suggested broadening mechanism 
experimentally. 

2. Doping process, e.g., random substitutions La — * 
Sr in La2- a; Sr z ,Cu04, intrinsically creates strong inho- 
mogencity in the system. By performing Hartrce-Fock 
calculations, we show that due to the strong Coulomb 
screening, ARPES lines obtain a very narrow width 
(r ~ 30 — 40 meV) as soon as doping is higher than 
~ 1%, in spite of the very strong disorder. These results 
provide a natural explanation for why the ARPES spec- 
tra undergo radical changes - from very broad Gaussian 
to narrow quasiparticle peaks - upon just a few percent 



doping of parent compounds. The residual "small" 
width r as 30 — 40 meV is due to the hole scattering from 
the residual (locally unscreened) part of the random 
potential. The residual width is quite universal, it is 
practically independent of doping/energy/momentum. 

3. The calculation of the surface density of states 
demonstrates that the top Cu02 layer of La2-xSr x Cu04 
is always in the Anderson localization regime, and we 
predict the Coulomb gap of the order of ~ 10 meV 
which could be observed with STM and/or ARPES 
experiments. 

4. The calculation of the bulk density of states also 
shows the Coulomb gap of the order of a few meV. 
The gap evolves smoothly through the percolation 
point x = 0.055. Hence the system remains in the 
Anderson localization regime, and this explains the 
insulating behavior observed in transport properties at 
high magnetic fields. 24 ' 25 

5. Considering Sr-doping induced disorder in 
La2- x Sr :Z ;Cu04, we find a very inhomogencous hole 
density profile which yields a broad NQR spectrum. 
The calculated doping and temperature dependencies of 
NQR lineshapes are consistent with experiments. 

Altogether, the results reported here highlight a signif- 
icant role played by Coulomb disorder effects in cuprates. 
In particular, screening of Coulomb defects (either of ex- 
trinsic origin or introduced by dopant ions) results in a 
dramatic evolution of physical properties upon doping. 
In this work, we focused mostly on the charge degrees 
of freedom, accounting for underlying magnetic correla- 
tions merely via a properly renormalized dispersion of the 
mobile holes. It remains a challenge to incorporate the 
magnetic degrees of freedom into the model explicitly, 
exploring thereby the coupled charge and spin dynamics 
in cuprates at short length scales. 
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